home *** CD-ROM | disk | FTP | other *** search
/ CD ROM Paradise Collection 4 / CD ROM Paradise Collection 4 1995 Nov.iso / hamradio / sgp4_pl2.zip / SGP_OBS.PAS < prev    next >
Pascal/Delphi Source File  |  1992-09-28  |  6KB  |  193 lines
  1. Unit SGP_Obs;
  2. {           Author:  Dr TS Kelso }
  3. { Original Version:  1992 Jun 02 }
  4. { Current Revision:  1992 Sep 28 }
  5. {          Version:  1.40 }
  6. {        Copyright:  1992, All Rights Reserved }
  7. {$N+}
  8.  
  9. INTERFACE
  10.   Uses SGP_Math;
  11.  
  12. Procedure Calculate_User_PosVel(var geodetic : vector;
  13.                                         time : double;
  14.                          var obs_pos,obs_vel : vector);
  15. Procedure Calculate_LatLonAlt(pos : vector;
  16.                              time : double;
  17.                      var geodetic : vector);
  18. Procedure Calculate_Obs(pos,vel,geodetic : vector;
  19.                                     time : double;
  20.                              var obs_set : vector);
  21. Procedure Calculate_RADec(pos,vel,geodetic : vector;
  22.                                       time : double;
  23.                                var obs_set : vector);
  24.  
  25. IMPLEMENTATION
  26.   Uses SGP_Intf,SGP_Init,SGP_Time;
  27.  
  28. Procedure Calculate_User_PosVel(var geodetic : vector;
  29.                                         time : double;
  30.                          var obs_pos,obs_vel : vector);
  31. { Reference:  The 1992 Astronomical Almanac, page K11. }
  32.   const
  33.     mfactor = twopi*omega_E/secday;
  34.   var
  35.     lat,lon,alt,
  36.     theta,c,s,achcp : double;
  37.   begin
  38.   lat := geodetic[1];
  39.   lon := geodetic[2];
  40.   alt := geodetic[3];
  41.   theta := Modulus(ThetaG_JD(time) + lon,twopi);
  42.   geodetic[4] := theta; {LMST}
  43.   c := 1/Sqrt(1 + f*(f - 2)*Sqr(Sin(lat)));
  44.   s := Sqr(1 - f)*c;
  45.   achcp := (xkmper*c + alt)*Cos(lat);
  46.   obs_pos[1] := achcp*Cos(theta);          {kilometers}
  47.   obs_pos[2] := achcp*Sin(theta);
  48.   obs_pos[3] := (xkmper*s + alt)*Sin(lat);
  49.   obs_vel[1] := -mfactor*obs_pos[2];       {kilometers/second}
  50.   obs_vel[2] :=  mfactor*obs_pos[1];
  51.   obs_vel[3] :=  0;
  52.   Magnitude(obs_pos);
  53.   Magnitude(obs_vel);
  54.   end; {Procedure Calculate_User_PosVel}
  55.  
  56. Procedure Calculate_LatLonAlt(pos : vector;
  57.                              time : double;
  58.                      var geodetic : vector);
  59. { Reference:  The 1992 Astronomical Almanac, page K12. }
  60.   var
  61.     lat,lon,alt,
  62.     theta,r,e2,phi,c : double;
  63.   begin
  64.   theta := AcTan(pos[2],pos[1]);
  65.   lon := Modulus(theta - ThetaG_JD(time),twopi);
  66.   r := Sqrt(Sqr(pos[1]) + Sqr(pos[2]));
  67.   e2 := f*(2 - f);
  68.   lat := AcTan(pos[3],r);
  69.   repeat
  70.     phi := lat;
  71.     c := 1/Sqrt(1 - e2*Sqr(Sin(phi)));
  72.     lat := AcTan(pos[3] + xkmper*c*e2*Sin(phi),r);
  73.   until Abs(lat - phi) < 1E-10;
  74.   alt := r/Cos(lat) - xkmper*c;
  75.   geodetic[1] := lat;   {radians}
  76.   geodetic[2] := lon;   {radians}
  77.   geodetic[3] := alt;   {kilometers}
  78.   geodetic[4] := theta; {radians}
  79.   end; {Procedure Calculate_LatLonAlt}
  80.  
  81. Procedure Calculate_Obs(pos,vel,geodetic : vector;
  82.                                     time : double;
  83.                              var obs_set : vector);
  84.   var
  85.     i                   : integer;
  86.     lat,lon,alt,theta,
  87.     sin_lat,cos_lat,
  88.     sin_theta,cos_theta : double;
  89.     el,azim             : double;
  90.     top_s,top_e,top_z   : double;
  91.     obs_pos,obs_vel,
  92.     range,rgvel         : vector;
  93.   begin
  94.   Calculate_User_PosVel(geodetic,time,obs_pos,obs_vel);
  95.   for i := 1 to 3 do
  96.     begin
  97.     range[i] := pos[i] - obs_pos[i];
  98.     rgvel[i] := vel[i] - obs_vel[i];
  99.     end; {for i}
  100.   Magnitude(range);
  101.   lat := geodetic[1];
  102.   lon := geodetic[2];
  103.   alt := geodetic[3];
  104.   theta := geodetic[4];
  105.   sin_lat := Sin(lat);
  106.   cos_lat := Cos(lat);
  107.   sin_theta := Sin(theta);
  108.   cos_theta := Cos(theta);
  109.   top_s := sin_lat*cos_theta*range[1]
  110.          + sin_lat*sin_theta*range[2]
  111.          - cos_lat*range[3];
  112.   top_e := -sin_theta*range[1]
  113.          + cos_theta*range[2];
  114.   top_z := cos_lat*cos_theta*range[1]
  115.          + cos_lat*sin_theta*range[2]
  116.          + sin_lat*range[3];
  117.   azim := ArcTan(-top_e/top_s); {Azimuth}
  118.   if top_s > 0 then
  119.     azim := azim + pi;
  120.   if azim < 0 then
  121.     azim := azim + twopi;
  122.   el := ArcSin(top_z/range[4]);
  123.   obs_set[1] := azim;                      {Azimuth (radians)}
  124.   obs_set[2] := el;                        {Elevation (radians)}
  125.   obs_set[3] := range[4];                  {Range (kilometers)}
  126.   obs_set[4] := Dot(range,rgvel)/range[4]; {Range Rate (kilometers/second)}
  127. { Corrections for atmospheric refraction }
  128. { Reference:  Astronomical Algorithms by Jean Meeus, pp. 101-104 }
  129. { Note:  Correction is meaningless when apparent elevation is below horizon }
  130.   obs_set[2] := obs_set[2] + Radians((1.02/Tan(Radians(Degrees(el)+10.3/(Degrees(el)+5.11))))/60);
  131.   if obs_set[2] >= 0 then
  132.     visible := true
  133.   else
  134.     begin
  135.     obs_set[2] := el;  {Reset to true elevation}
  136.     visible := false;
  137.     end; {else}
  138.   end; {Procedure Calculate_Obs}
  139.  
  140. Procedure Calculate_RADec(pos,vel,geodetic : vector;
  141.                                       time : double;
  142.                                var obs_set : vector);
  143. { Reference:  Methods of Orbit Determination by Pedro Ramon Escobal, pp. 401-402}
  144.   var
  145.     phi,theta,
  146.     sin_theta,cos_theta,
  147.     sin_phi,cos_phi,
  148.     az,el,
  149.     Lxh,Lyh,Lzh,
  150.     Sx,Ex,Zx,
  151.     Sy,Ey,Zy,
  152.     Sz,Ez,ZZ,
  153.     Lx,Ly,Lz,
  154.     cos_delta,
  155.     sin_alpha,cos_alpha : double;
  156.   begin
  157.   Calculate_Obs(pos,vel,geodetic,time,obs_set);
  158.   if visible then
  159.     begin
  160.     az := obs_set[1];
  161.     el := obs_set[2];
  162.     phi   := geodetic[1];
  163.     theta := Modulus(ThetaG_JD(time) + geodetic[2],twopi);
  164.     sin_theta := Sin(theta);
  165.     cos_theta := Cos(theta);
  166.     sin_phi := Sin(phi);
  167.     cos_phi := Cos(phi);
  168.     Lxh := -Cos(az)*Cos(el);
  169.     Lyh :=  Sin(az)*Cos(el);
  170.     Lzh :=  Sin(el);
  171.     Sx := sin_phi*cos_theta;
  172.     Ex := -sin_theta;
  173.     Zx := cos_theta*cos_phi;
  174.     Sy := sin_phi*sin_theta;
  175.     Ey := cos_theta;
  176.     Zy := sin_theta*cos_phi;
  177.     Sz := -cos_phi;
  178.     Ez := 0;
  179.     Zz := sin_phi;
  180.     Lx := Sx*Lxh + Ex*Lyh + Zx*Lzh;
  181.     Ly := Sy*Lxh + Ey*Lyh + Zy*Lzh;
  182.     Lz := Sz*Lxh + Ez*Lyh + Zz*Lzh;
  183.     obs_set[2] := ArcSin(Lz);                 {Declination (radians)}
  184.     cos_delta := Sqrt(1 - Sqr(Lz));
  185.     sin_alpha := Ly/cos_delta;
  186.     cos_alpha := Lx/cos_delta;
  187.     obs_set[1] := AcTan(sin_alpha,cos_alpha); {Right Ascension (radians)}
  188.     obs_set[1] := Modulus(obs_set[1],twopi);
  189.     end; {if}
  190.   end; {Procedure Calculate_RADec}
  191.  
  192. end.
  193.